Using R as a Geographic Information System

Toyin Ola

Why use R as a GIS?

Example: Finding Maternity Care Deserts

Consider what steps you would have to take if you wanted to see which parts of AZ do not have a hospital with inpatient labor and delivery (L&D) units within a 45-minute drive.

  1. Obtain a list of hospitals with L&D.
  2. Geocode hospital addresses.
  3. Determine a 45-minute radius around each hospital.
  4. Obtain AZ geographic boundaries.
  5. Visualize the results.

Set Up

Package Purpose
jsonlite Query CMS API to obtain data on hospitals
tidygeocoder Geocode hospital addresses
dplyr Wrangle data
stringr Use text to filter data
tigris Download Census Bureau TIGER/Line geographic boundaries
sf Manipulate geosptial data
hereR Query HERE API to determine 45-minute driving radius
mapview Visualize geospatial data quickly
leaflet Generate interactive geospatial visualizations

To help query the HERE API, there is one user-defined tryLocation() function to source. This function is courtesy of Hannah Recht.


Step 1

Obtain a list of AZ hospitals with L&D from the CMS Provider Data Catalog. Determine which hospitals do not have L&D.

# Query CMS Provider Data Catalog

query <- paste0(
  "https://data.cms.gov/provider-data/api/1/datastore/query/nrdb-3fcy/0",
  "?conditions[0][property]=State&conditions[0][value]=", params$state,
  "&conditions[0][operator]=="
)
  
raw <- jsonlite::fromJSON(query)
  
hospital_df <- raw$results

# Identify hospitals w/out L&D

no_inpatient <- hospital_df |>  
  filter(stringr::str_detect(score, "does not provide inpatient labor")) |> 
  distinct(facility_id, .keep_all = TRUE) |> 
  select(c(facility_id, facility_name)) 

# Preview 

head(hospital_df)

head(no_inpatient)
  facility_id                              facility_name
1       30002 BANNER - UNIVERSITY MEDICAL CENTER PHOENIX
2       30002 BANNER - UNIVERSITY MEDICAL CENTER PHOENIX
3       30002 BANNER - UNIVERSITY MEDICAL CENTER PHOENIX
4       30002 BANNER - UNIVERSITY MEDICAL CENTER PHOENIX
5       30002 BANNER - UNIVERSITY MEDICAL CENTER PHOENIX
6       30002 BANNER - UNIVERSITY MEDICAL CENTER PHOENIX
                  address citytown state zip_code countyparish telephone_number
1 1111 EAST MCDOWELL ROAD  PHOENIX    AZ    85006     MARICOPA   (602) 839-2000
2 1111 EAST MCDOWELL ROAD  PHOENIX    AZ    85006     MARICOPA   (602) 839-2000
3 1111 EAST MCDOWELL ROAD  PHOENIX    AZ    85006     MARICOPA   (602) 839-2000
4 1111 EAST MCDOWELL ROAD  PHOENIX    AZ    85006     MARICOPA   (602) 839-2000
5 1111 EAST MCDOWELL ROAD  PHOENIX    AZ    85006     MARICOPA   (602) 839-2000
6 1111 EAST MCDOWELL ROAD  PHOENIX    AZ    85006     MARICOPA   (602) 839-2000
  measure_id
1      PC_01
2      PC_05
3       SM_7
4     ePC_02
5    ePC_07a
6    ePC_07b
                                                                           measure_name
1                                                                     Elective Delivery
2                                                         Exclusive Breast Milk Feeding
3                                                 Maternal Morbidity Structural Measure
4                                                                        Cesarean Birth
5                                    Risk Adjusted Severe Obstetric Complications (All)
6 Risk Adjusted Severe Obstetric Complications (excluding blood-transfusion-only cases)
          score         sample footnote start_date   end_date
1             1            109        2 01/01/2023 12/31/2023
2 Not Available  Not Available        5 01/01/2023 12/31/2023
3           Yes Not Applicable          01/01/2023 12/31/2023
4 Not Available  Not Available        5 01/01/2023 12/31/2023
5 Not Available  Not Available        5 01/01/2023 12/31/2023
6 Not Available  Not Available        5 01/01/2023 12/31/2023
  facility_id                                facility_name
1       30010                          ST. MARY'S HOSPITAL
2       30014  HONOR HEALTH JOHN C. LINCOLN MEDICAL CENTER
3       30030                        ABRAZO CENTRAL CAMPUS
4       30037                      TEMPE ST LUKES HOSPITAL
5       30038 HONORHEALTH SCOTTSDALE OSBORN MEDICAL CENTER
6       30061                BANNER BOSWELL MEDICAL CENTER

Step 2

Geocode the hospitals addresses and turn into a geospatial dataset.

# Filter and reformat hospital data

hospital_df <- hospital_df |> 
  distinct(facility_id, .keep_all = TRUE) |> 
  mutate(inpatient = case_when(facility_id %in% no_inpatient$facility_id  ~ "No",
                               .default = "Yes")) |> 
  mutate(full_address = paste0(address, ", ", citytown, ", ", state, " ", zip_code),
         .after = zip_code) |> 
  select(-c(telephone_number:end_date))

# Geocode hospital addresses

hospital_df <- hospital_df |> 
  tidygeocoder::geocode(address = full_address, method = "arcgis") 

# Convert to sf object

hospitals <- hospital_df |>
  st_as_sf(coords = c("long", "lat"), 
           crs = 4326) # use {leaflet} preferred CRS

# Preview for quick QC

mapview::mapview(hospitals)

Step 3

Query the HERE API to determine a 45-minute driving radius (i.e., isochrone) around each hospital.

# Set HERE API key

set_key(Sys.getenv("HERE_API_KEY"))

# Loop over point data to make isochrones file (using Hannah Recht's function)

## select hospitals w/ L&D

ip_hospitals <- hospitals |> 
  filter(inpatient == "Yes")

## create empty vectors for loop output

isochrones <- vector(mode = "list", length = nrow(ip_hospitals))
error_rows <- vector(mode = "list", length = nrow(ip_hospitals))

for (i in 1:nrow(ip_hospitals)) {
  
    print(i)
  
    ## get isochrone for that point, using delay to avoid rate limiting
    Sys.sleep(0.9)
    
    ## filter to ith point
    point_temp <- ip_hospitals %>% filter(row_number() == i)
    point_id <- point_temp$facility_id
    
    isochrones_temp <- tryLocation(point_temp)
    
    ## save any errored out points
    if (is.null(isochrones_temp)) {
        error_rows <- bind_rows(error_rows, point_temp)
    } else {
        isochrones <- bind_rows(isochrones, isochrones_temp)    
    }
}

# Remove extraneous columns

isochrones <- isochrones |> 
  select(facility_id, geometry)

# Left join hospital data

ip_hospitals <- ip_hospitals |> 
  select(c(facility_id, facility_name, full_address, inpatient)) |> 
  st_drop_geometry()

isochrones <- isochrones |> 
  left_join(ip_hospitals, by = join_by(facility_id))

# Reproject

isochrones <- isochrones |> 
  st_transform(crs = 4326)

# Preview

mapview::mapview(isochrones)

If you want to try it out without signing up for a HERE API key, use the data available on GitHub.

Step 4

Obtain AZ county boundaries.

# Download county boundaries 

counties <- counties(state = c(params$state), 
                     cb = TRUE, 
                     resolution = "500k", 
                     year = 2023) 

# Reproject to {leaflet} preferred CRS

counties <- counties |> 
  st_transform(crs = 4326)

# Create a quick QC viz

mapview::mapview(counties)

Step 5

Visualize the hospitals and isochrones in an interactive map.

# Specify text for popups

hosp_popup <- paste0(
  "<b>", hospitals$facility_name, "</b>", "<br/>",
  hospitals$full_address
)

iso_popup <- paste0(
  "This is a 45-minute driving radius around ", isochrones$facility_name, "."
)

# Create color palette

color_icon <- awesomeIcons(
  icon = 'location-dot',
  iconColor = 'black',
  library = 'fa',
  markerColor = ifelse(hospitals$inpatient == 'Yes', 'green', 'gray')) # no HEX codes

# Make an interactive map

leaflet() %>%
  addProviderTiles("CartoDB.Voyager") %>%    # add base map
  addPolygons(data = counties,
              weight = 1,
              color = "#999999",
              fillColor = "#d4d4d4",
              stroke = TRUE,
              popup = paste0(counties$NAME)) %>%
  addAwesomeMarkers(data = hospitals, 
             icon = color_icon,
             group = hospitals$facility_name,
             popup = ~hosp_popup) %>%
  addPolygons(data = isochrones,
              weight = 1, opacity = 1.0,     # set stroke width and opacity
              stroke = TRUE,
              color = "#669933",
              fillColor = "#669933", 
              fillOpacity = 0.5, 
              group = isochrones$facility_name,
              popup = ~iso_popup,
              highlightOptions = highlightOptions(color = "#527c29",
                                                  weight = 3,
                                                  bringToFront = TRUE)
  ) %>%
  addLayersControl(
    overlayGroups = sort(c(isochrones$facility_name, hospitals$facility_name)),
    options = layersControlOptions(collapsed = TRUE)
  ) 

Possibilities

Learn More

Here are some resources I have found helpful translating knowledge from a GIS like QGIS to R: